##########################################################################
#             REPLICATION CODE FOR:                                      #
# Foreign Sponsors and Inter-Rebel Alliances                             #
# Milos Popovic, 11/9/2017                                               #
#                                                                        #
##########################################################################
# R version 3.4.1 (2017-06-30) -- "Single Candle"						 #
##########################################################################
### BEGIN REPLICATION													 #
##########################################################################

install.packages("gridExtra",dependencies=TRUE)
install.packages("ggplot2",dependencies=TRUE)
install.packages("R2admb")
install.packages("glmmADMB", 
       repos=c("http://glmmadmb.r-forge.r-project.org/repos",
       getOption("repos")),
       type="source")
install.packages("Amelia",dependencies=TRUE)
install.packages("mi",dependencies=TRUE)
install.packages("BBmisc",dependencies=TRUE)
install.packages("Hmisc",dependencies=TRUE)
install.packages("plyr",dependencies=TRUE)
install.packages("data.table",dependencies=TRUE)
install.packages("ggthemes",dependencies=TRUE)
install.packages("plyr",dependencies=TRUE)
install.packages("dplyr",dependencies=TRUE)
install.packages("ROCR",dependencies=TRUE)
install.packages("SDMTools",dependencies=TRUE)
install.packages("mitools",dependencies=TRUE)
install.packages("devtools",dependencies=TRUE)
install.packages("ModelCriticism",dependencies=TRUE)
install.packages("sampleSelection",dependencies=TRUE)
install.packages("cvTools",dependencies=TRUE)
install.packages("caret",dependencies=TRUE)
install.packages("parallel",dependencies=TRUE)
install_github("zsmahmood89/ModelCriticism/packages/ModelCriticism")

library(grid)
library(gridExtra)
library(ggplot2)
library(glmmADMB)
library(R2admb)
library(Amelia)
library(mi)
library(BBmisc)
library(Hmisc)
library(plyr)
library(data.table)
library(ggthemes)
library(dplyr)
library(ROCR)
library(SDMTools)
library(mitools)
library(devtools)
library(ModelCriticism)
library(sampleSelection)
library(cvTools)
library(caret)
library(parallel)

# prepare data

a <- read.csv("final.csv", header=T)
var_out <- c("insurgentsa", "insurgentsb", 'size_a',	'size_b',	'id',	'location',	'oldid',	'target', 'dur', 'statesupport')

a$dyad <- paste(a$insurgentsa,a$insurgentsb, sep = ' & ') # create dyad labels
a$dyadid <- as.character(as.numeric(factor(a$dyad,levels=unique(a$dyad)))) # create dyad numbers
#print(a$dyad)
a$dyad <- toupper(a$dyad) # capitalize dyad labels

d <- a[, !(names(a)%in%var_out)] 

#frequency plot or FIGURE 1 IN THE ARTICLE

nrow(a[a$sponsor==0 & a$ally==1,])
nrow(a[a$sponsor==1 & a$ally==1,])
nrow(a[a$sponsor==2 & a$ally==1,])
nrow(a[a$sponsor==3 & a$ally==1,])
nrow(a[a$sponsor==0 & a$ally==0,])
nrow(a[a$sponsor==1 & a$ally==0,])
nrow(a[a$sponsor==2 & a$ally==0,])
nrow(a[a$sponsor==3 & a$ally==0,])

tab = as.table(rbind(c(90, 26, 5, 85),
                     c(607, 96, 14, 62)))
names(dimnames(tab)) = c("Alliance", "Sponsor")
rownames(tab)        = c("Yes", "No")
colnames(tab)        = c("None", "Single", "Different", "Shared")
colors <- c("grey80", "grey30")

gr <- (1+sqrt(5))/2 # golden ration for plots

pdf('fig1.pdf', 10, 10*(gr^-1))
par(xpd=F, mar=c(5,5,1,1))
mp <- barplot(tab, beside = TRUE, main="", xlab = "Sponsor", ylab = "Number of observations", ylim=c(0,max(tab)+50), col=colors, args.legend = list(x="center", fill=colors, title = "Alliance?", x=ncol(tab) + 3, y=max(colSums(tab))))
text(mp, tab, labels = format(tab, 4),
pos = 3, cex = .75)
legend(x=4, y=600, rownames(tab), fill = colors, bty = "n", title="Alliance?", xpd=NA)
dev.off()

# imputation

apply(d, 2, function(i) mean(is.na(i)))

d$cowcode <- as.factor(d$cowcode)

d$ratio[d$ratio==1] <- 1-1e-3
d$ratio[d$ratio==0] <- 0+1e-3

d$relfrac[d$relfrac==1] <- 1-1e-3
d$relfrac[d$relfrac==0] <- 0+1e-3

d$ef[d$ef==1] <- 1-1e-3
d$ef[d$ef==0] <- 0+1e-3

set.seed(1914)

b <- amelia(x=d, m=5, 
ts='year', 
cs='cowcode', 
lgstc='ratio', 
noms='wd', 
idvars=c('dyad', 'dyadid'),
bounds=rbind(c(8, 0, 1), c(9, 0, 1), c(10, 0, Inf), c(11, 0, Inf), c(16, 0, Inf)),
polytime=3)

#### prepare models####

fitb <- ally~ 1+lgdppc+lnmilex+ldur+ncontig+relfrac+ef+lndurable+(1|cowcode) #base fit not shown in the paper
fit0 <- ally~ 1+rsup*wd+lgdppc+lnmilex+ldur+ncontig+relfrac+ef+lndurable+(1|cowcode)
fit0b <- ally~ 1+ratio+lgdppc+lnmilex+ldur+ncontig+relfrac+ef+lndurable+(1|cowcode)
fit <- ally~ 1+factor(sponsor)+lgdppc+lnmilex+ldur+ncontig+relfrac+ef+lndurable+(1|cowcode)

### function for running glmmADMB on MIs
burrito <- function(x, form){ glmmadmb(form, data=x, family='binomial')
}

### fit models 

g0 <- lapply(b$imputations, burrito, form=fit0)
g0b <- lapply(b$imputations, burrito, form=fit0b)
g <- lapply(b$imputations, burrito, form=fit)
base <- lapply(b$imputations, burrito, form=fitb) #base fit not shown in the paper

### Save for future use

save(list=c('g', 'g0', 'g0b', 'base'), file='models.RData')

# pool coefficients for model 1 by Rubin's rule

gcoefs <- do.call(rbind, lapply(g, function(i) coef(summary(i))[,1])) # coef estimates
gses <- do.call(rbind, lapply(g, function(i) coef(summary(i))[,2])) # standard errors
gd <- data.frame(t(do.call(rbind, mi.meld(gcoefs, gses)))) # pool coefs and errors by Rubin's rule 
names(gd) <- c("coef", "se")
gd <- gd[-c(1), ] # get rid of intercept
gd$lo <- gd$coef - 1.96*gd$se # lower 95% confidence intervals
gd$up <- gd$coef + 1.96*gd$se # upper 95% confidence intervals
gd$var = c("Single sponsor", "Different sponsor", "Shared sponsor", "GDP p.c. (ln)", "Expenditure (ln)", "Duration (ln)", "Non-Contiguity", "Religious frac.", "Ethnic frac.", "Durability (ln)") # assign names to vars

# plot model 1

p <- ggplot(data=gd) + scale_x_discrete(limits = c("Durability (ln)", "Ethnic frac.", "Religious frac.", "Non-Contiguity", "Duration (ln)", "Expenditure (ln)", "GDP p.c. (ln)", "Shared sponsor", "Different sponsor",   "Single sponsor"))
p <- p + geom_hline(yintercept = 0, colour = "grey20", lty = 2) 
p <- p + geom_linerange(aes(x = var, ymin = lo,
                            ymax = up), col= "grey20",
                            lwd = 1, position = position_dodge(width = 1/2))
p <- p + geom_pointrange(aes(x = var, y = coef, ymin = lo,
                                 ymax = up), col= "grey20",
                             lwd = 1/2, position = position_dodge(width = 1/2),
                             shape = 21, fill = "grey20")
p <- p + coord_flip() + guides(fill=FALSE) + theme(panel.background = element_blank(),
              panel.grid.minor = element_blank(),
              panel.grid.major=element_blank(),
              panel.border=element_blank(),
			  plot.title=element_text(hjust = .5, size=11),
              axis.title.x = element_text(colour="grey30"),
              axis.title.y = element_text(colour="grey30"),
              axis.line   = element_line(colour=NA),
              axis.line.x = element_line(colour="grey30"),
              axis.line.y = element_line(colour="grey30"),
              axis.text.x = element_text(colour="grey30", hjust = 1),
              axis.text.y = element_text(colour="grey30"), legend.title = element_blank(),
              legend.key=element_rect(fill=NA),
              legend.position="bottom") + ylab("Coefficient Estimates") + xlab(NULL) +
			  ggtitle("Model 1: Alliance in the Shadow\n of Sponsor")

# pool coefficients for model 2 by Rubin's rule

gcoefs0 <- do.call(rbind, lapply(g0, function(i) coef(summary(i))[,1]))
gses0 <- do.call(rbind, lapply(g0, function(i) coef(summary(i))[,2]))
gd1 <- data.frame(t(do.call(rbind, mi.meld(gcoefs0, gses0))))
names(gd1) <- c("coef", "se") #change names
gd1 <- gd1[-c(1), ] # get rid of intercept
gd1$lo <- gd1$coef - 1.96*gd1$se # lower 95% confidence intervals
gd1$up <- gd1$coef + 1.96*gd1$se # upper 95% confidence intervals
gd1$var = c("Sponsor", "Weak Dyad", "GDP p.c. (ln)", "Expenditure (ln)", "Duration (ln)", "Non-Contiguity", "Religious frac.", "Ethnic frac.", "Durability (ln)", "Sponsor x\n Weak Dyad")

# plot model 2

p2 <- ggplot(data=gd1) + scale_x_discrete(limits = c("Durability (ln)", "Ethnic frac.", "Religious frac.", "Non-Contiguity", "Duration (ln)", "Expenditure (ln)", "GDP p.c. (ln)","Sponsor x\n Weak Dyad", "Sponsor", "Weak Dyad"))
p2 <- p2 + geom_hline(yintercept = 0, colour = "grey20", lty = 2) 
p2 <- p2 + geom_linerange(aes(x = var, ymin = lo,
                            ymax = up), col= "grey20",
                            lwd = 1, position = position_dodge(width = 1/2))
p2 <- p2 + geom_pointrange(aes(x = var, y = coef, ymin = lo,
                                 ymax = up), col= "grey20",
                             lwd = 1/2, position = position_dodge(width = 1/2),
                             shape = 21, fill = "grey20")
p2 <- p2 + coord_flip() + guides(fill=FALSE) + theme(panel.background = element_blank(),
              panel.grid.minor = element_blank(),
              panel.grid.major=element_blank(),
              panel.border=element_blank(),
			  plot.title=element_text(hjust = .5, size=11),
              axis.title.x = element_text(colour="grey30"),
              axis.title.y = element_text(colour="grey30"),
              axis.line   = element_line(colour=NA),
              axis.line.x = element_line(colour="grey30"),
              axis.line.y = element_line(colour="grey30"),
              axis.text.x = element_text(colour="grey30", hjust = 1),
              axis.text.y = element_text(colour="grey30"), legend.title = element_blank(),
              legend.key=element_rect(fill=NA),
              legend.position="bottom") + ylab("Coefficient Estimates") + xlab(NULL) +
			  ggtitle("Model 2: Sponsor and Weak Dyad\n (Bapat and Bond)")

# pool coefficients for model 3 by Rubin's rule

gcoefs0b <- do.call(rbind, lapply(g0b, function(i) coef(summary(i))[,1]))
gses0b <- do.call(rbind, lapply(g0b, function(i) coef(summary(i))[,2]))
gd2 <- data.frame(t(do.call(rbind, mi.meld(gcoefs0b, gses0b))))
names(gd2) <- c("coef", "se")
gd2 <- gd2[-c(1), ] # get rid of intercept
gd2$lo <- gd2$coef - 1.96*gd2$se # lower 95% confidence intervals
gd2$up <- gd2$coef + 1.96*gd2$se # upper 95% confidence intervals
gd2$var = c("Ratio", "GDP p.c. (ln)", "Expenditure (ln)", "Duration (ln)", "Non-Contiguity", "Religious frac.", "Ethnic frac.", "Durability (ln)")

# plot model 3

p3 <- ggplot(data=gd2) + scale_x_discrete(limits = c("Durability (ln)", "Ethnic frac.", "Religious frac.", "Non-Contiguity", "Duration (ln)", "Expenditure (ln)", "GDP p.c. (ln)", "Ratio"))
p3 <- p3 + geom_hline(yintercept = 0, colour = "grey20", lty = 2) 
p3 <- p3 + geom_linerange(aes(x = var, ymin = lo,
                            ymax = up), col= "grey20",
                            lwd = 1, position = position_dodge(width = 1/2))
p3 <- p3 + geom_pointrange(aes(x = var, y = coef, ymin = lo,
                                 ymax = up), col= "grey20",
                             lwd = 1/2, position = position_dodge(width = 1/2),
                             shape = 21, fill = "grey20")
p3 <- p3 + coord_flip() + guides(fill=FALSE) + theme(panel.background = element_blank(),
              panel.grid.minor = element_blank(),
              panel.grid.major=element_blank(),
              panel.border=element_blank(),
			  plot.title=element_text(hjust = .5, size=11),
              axis.title.x = element_text(colour="grey30"),
              axis.title.y = element_text(colour="grey30"),
              axis.line   = element_line(colour=NA),
              axis.line.x = element_line(colour="grey30"),
              axis.line.y = element_line(colour="grey30"),
              axis.text.x = element_text(colour="grey30", hjust = 1),
              axis.text.y = element_text(colour="grey30"), legend.title = element_blank(),
              legend.key=element_rect(fill=NA),
              legend.position="bottom") + ylab("Coefficient Estimates") + xlab(NULL) +
			  ggtitle("Model 3: Relative Strength (Christia)")

# COMBINE MODELS INTO A SINGLE PLOT (FIGURE 2 IN THE ARTICLE)
			  
pp <- grid.arrange(p, p2, p3, ncol=2, nrow=2)			  
ggsave("fig2.pdf", height=7, width=7, pp)

# FIGURE 3 IN THE ARTICLE
#predicted probabilites as "Sponsor" varies and other predictors are held constant
#prepare new data frame composed of mean or median values

nd <- expand.grid(ally=factor(0:1),
sponsor=factor(rep(0:3, each = 2)), 
lgdppc = 2.5,
ncontig = 0,
lnmilex = 12,
ldur = 1.14,
relfrac = 0.48,
ef = 0.71,
lndurable = 1.83)

#pool predictions with standard errors
p_eta <- lapply(g, predict, newdata=nd, re.form=NA, type = "link", se = TRUE)
pr_eta <- do.call(rbind, p_eta)
l = apply( cbind( pr_eta ) , 1 , unlist )
fdd <- t(l[c(1:16),]) # probabilities
sdd <- t(l[c(17:32),] ) # standard errors
p_invlog <- plogis(fdd)[ ,c(1:16)] # log-transform
quants_invlog <- apply(p_invlog, 2, quantile, probs = c(0, .25, .5, .75, 1)) #generate quantiles for boxplot
df <- cbind(nd, t(quants_invlog))
df1 <- ddply(df, .(sponsor), summarise, ymin = mean(`0%`), lower = mean(`25%`), middle = mean(`50%`), upper = mean(`75%`), ymax = mean(`100%`))

# plot
p1 <- ggplot(df1, aes(x=sponsor, y=middle)) +
  geom_boxplot(
   aes(ymin = ymin, lower = lower, middle = middle, upper = upper, ymax = ymax),
   stat = "identity"
 ) + scale_x_discrete(labels = c("No Sponsor", "Single Sponsor", "Different Sponsors",  "Shared Sponsor")) + 
scale_colour_manual(values=c("grey70", "grey50", "grey30", 
                                 "grey10")) +
                    guides(fill=FALSE) + 
theme(panel.background = element_blank(),
              panel.grid.minor = element_blank(),
              panel.grid.major=element_blank(),
              panel.border=element_blank(),
              axis.title.x = element_text(colour="grey30"),
              axis.title.y = element_text(colour="grey30"),
              axis.line   = element_line(colour=NA),
              axis.line.x = element_line(colour="grey30"),
              axis.line.y = element_line(colour="grey30"),
              axis.text.x = element_text(colour="grey30", hjust = 0.5),
              axis.text.y = element_text(colour="grey30"), legend.title = element_blank(),
              legend.key=element_rect(fill=NA),
              legend.position="none") + ylab("Pr Alliance\n") + xlab("Sponsor\n")

ggsave("fig3.pdf", height=7, width=7, p1)

# predicted probabilities Christia
nd2 <- expand.grid(ally=factor(0:1),
ratio=rep(seq(from=0, to=1, length.out = 100)), 
lgdppc = 2.5,
ncontig = 0,
lnmilex = 12,
ldur = 1.14,
relfrac = 0.48,
ef = 0.71,
lndurable = 1.77)

p_eta2 <- lapply(g0b, predict, newdata=nd2, re.form=NA, type = "link", se = TRUE)
pr_eta2 <- do.call(rbind, p_eta2)
l = apply( cbind( pr_eta2 ) , 1 , unlist )
fdd <- t(l[c(1:200),]) 
sdd <- t(l[c(201:400),] )
da <- data.frame(t(do.call(rbind, mi.meld(fdd, sdd))))
names(da)[1] <- "prob"
names(da)[2] <- "se"
da$left <- plogis(da$prob-(1.96*da$se))
da$right <- plogis(da$prob+(1.96*da$se))
da$fit <- plogis(da$prob)
dfd <- data.frame(da, nd2)
  
q1 <- ggplot(dfd, aes(x = ratio, y = fit)) +  
geom_line(size = 1) + 
geom_ribbon(aes(ymin=left, ymax=right), alpha=0.2) +
theme_bw() +
theme(panel.background = element_blank(),
              panel.grid.minor = element_blank(),
              panel.grid.major=element_blank(),
              panel.border=element_blank(),
              axis.title.x = element_text(colour="grey30"),
              axis.title.y = element_text(colour="grey30"),
              axis.line   = element_line(colour=NA),
              axis.line.x = element_line(colour="grey30"),
              axis.line.y = element_line(colour="grey30"),
              axis.text.x = element_text(colour="grey30", hjust = 0.5),
              axis.text.y = element_text(colour="grey30"), legend.title = element_blank(),
              legend.key=element_rect(fill=NA),
              legend.position="none") + ylab("") + xlab("Ratio")
  
  ppp <- grid.arrange(p1,q1,nrow=1, ncol=2)
  ggsave("fig4.pdf", height=7, width=10.5, ppp)

# 10-FOLD CROSS-VALIDATION

#Make a new dataframe of necessary variables to put values into
f <- d[,c("ally", "cowcode", "sponsor", 
"rsup", "wd", "ratio", "lgdppc", "lnmilex", "ldur", "ncontig", "relfrac","ef", "lndurable")]

f$cowcode <- as.factor(f$cowcode)
#f$sponsor <- as.factor(f$sponsor)

# prediction and testset data frames that we add to with each iteration over
# the 10 folds
# prepare the models

fit1 <- ally~ 1+factor(sponsor)+lgdppc+lnmilex+ldur+ncontig+relfrac+ef+lndurable+(1|cowcode)

fit2 <- ally~ 1+rsup*wd+lgdppc+lnmilex+ldur+ncontig+relfrac+ef+lndurable+(1|cowcode)

fit3 <- ally~ 1+ratio+lgdppc+lnmilex+ldur+ncontig+relfrac+ef+lndurable+(1|cowcode)

# sample from 1 to k, nrow times (the number of observations in the data)
#   this needs a seed so it's reproducible
#   the folds are of unequal sizes
#K-fold cross-validation
set.seed(667)
n_folds <- 1e1 
nrow(f)/n_folds
for_folds <- rep(1:n_folds, c(99, 98, 99, 98, 99,  98, 99, 98, 99, 98))
nrow(f)==length(for_folds)
f$id <- sample(for_folds, nrow(f), replace = FALSE)
id_vec <- 1:n_folds

# Do multiple imputation

#   Amelia has trouble with the logistic for f$ratio, bc 17 values equal 1
#   here's the conventional hack of making it very slightly < 1 
f$ratio[f$ratio==1] <- 1-1e-12

set.seed(1914)
v <- amelia(x=f, m=5, cs='cowcode', lgstc='ratio', noms='wd', idvars='sponsor',
            bounds=rbind(c(5, 0, 1), c(6, 0, 1), c(7, 0, Inf), c(8, 0, Inf), c(13, 0, Inf))) 


set.seed(1989)
#  the filter for the train set needs to be in!
pizza <- 
    function(x, form, fold_id) 
    { 
        #   trainf <- subset(x, id %in% id_vec[-fold_id])
        trainf <- x[f$id!=fold_id, ]
        #   testf <- subset(f, id %in% c(i))
        glmmadmb(form, data=trainf, family='binomial') 
    }

# function to apply CV on multiply imputed data

pizzamod1 <-
    function(i) 
    {
        lapply(v$imputations, pizza, form=fit1, fold_id=i) 
    }
	
pizzamod2 <-
    function(i) 
    {
        lapply(v$imputations, pizza, form=fit2, fold_id=i)
    }
	
pizzamod3 <-
    function(i) 
    {
        lapply(v$imputations, pizza, form=fit3, fold_id=i)
    }
	
system.time(list_m1 <- lapply(id_vec, pizzamod1)) # imputations for Model 1
system.time(list_m2 <- lapply(id_vec, pizzamod2)) # imputations for Model 2
system.time(list_m3 <- lapply(id_vec, pizzamod3)) # imputations for Model 3

save(list=c('f', 'v', 'list_m1', 'list_m2', 'list_m3'), file='cv2.RData')

# extract data for each fold from 500 imputations

ids <- lapply(1:5, function(i) subset(v$imputations[[i]], f$id ==1))
ids1 <- lapply(1:5, function(i) subset(v$imputations[[i]], f$id ==2))
ids2 <- lapply(1:5, function(i) subset(v$imputations[[i]], f$id ==3))
ids3 <- lapply(1:5, function(i) subset(v$imputations[[i]], f$id == 4))
ids4 <- lapply(1:5, function(i) subset(v$imputations[[i]], f$id ==5))
ids5 <- lapply(1:5, function(i) subset(v$imputations[[i]], f$id == 6))
ids6 <- lapply(1:5, function(i) subset(v$imputations[[i]], f$id ==7))
ids7 <- lapply(1:5, function(i) subset(v$imputations[[i]], f$id == 8))
ids8 <- lapply(1:5, function(i) subset(v$imputations[[i]], f$id ==9))
ids9 <- lapply(1:5, function(i) subset(v$imputations[[i]], f$id == 10))

# MODEL 1

coefs <- do.call(rbind, lapply(list_m1[[1]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m1[[1]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids)
newX <- model.matrix(list_m1[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
mu <- as.data.frame(plogis(eta))
names(mu)[1] <- "fit"
mu$id <- 1

coefs <- do.call(rbind, lapply(list_m1[[2]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m1[[2]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids1)
newX <- model.matrix(list_m1[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
mu1 <- as.data.frame(plogis(eta))
names(mu1)[1] <- "fit"
mu1$id <- 2

coefs <- do.call(rbind, lapply(list_m1[[3]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m1[[3]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids2)
newX <- model.matrix(list_m1[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
mu2 <- as.data.frame(plogis(eta))
names(mu2)[1] <- "fit"
mu2$id <- 3

coefs <- do.call(rbind, lapply(list_m1[[4]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m1[[4]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
d1 <- d1[-3] #else: returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids3)
newX <- model.matrix(list_m1[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
mu3 <- as.data.frame(plogis(eta))
names(mu3)[1] <- "fit"
mu3$id <- 4

coefs <- do.call(rbind, lapply(list_m1[[5]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m1[[5]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids4)
newX <- model.matrix(list_m1[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
mu4 <- as.data.frame(plogis(eta))
names(mu4)[1] <- "fit"
mu4$id <- 5

coefs <- do.call(rbind, lapply(list_m1[[6]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m1[[6]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids5)
newX <- model.matrix(list_m1[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
mu5 <- as.data.frame(plogis(eta))
names(mu5)[1] <- "fit"
mu5$id <- 6

coefs <- do.call(rbind, lapply(list_m1[[7]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m1[[7]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
d1 <- as.vector(d1)
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
nd <- as.data.frame(ids6)
newX <- model.matrix(list_m1[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
mu6 <- as.data.frame(plogis(eta))
names(mu6)[1] <- "fit"
mu6$id <- 7

coefs <- do.call(rbind, lapply(list_m1[[8]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m1[[8]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids7)
newX <- model.matrix(list_m1[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
mu7 <- as.data.frame(plogis(eta))
names(mu7)[1] <- "fit"
mu7$id <- 8

coefs <- do.call(rbind, lapply(list_m1[[9]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m1[[9]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids8)
newX <- model.matrix(list_m1[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
mu8 <- as.data.frame(plogis(eta))
names(mu8)[1] <- "fit"
mu8$id <- 9

coefs <- do.call(rbind, lapply(list_m1[[10]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m1[[10]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids9)
newX <- model.matrix(list_m1[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
mu9 <- as.data.frame(plogis(eta))
names(mu9)[1] <- "fit"
mu9$id <- 10

 list_yhat <- list(mu, mu1, mu2, mu3, mu4, mu5,mu6,mu7,mu8,mu9)
 
 
for (i in 1:10) {
f$yhat[f$id==i] <- list_yhat[[i]]$fit
}

# MODEL 2

coefs <- do.call(rbind, lapply(list_m2[[1]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m2[[1]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids)
newX <- model.matrix(list_m2[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pr <- as.data.frame(plogis(eta))
names(pr)[1] <- "fit"
pr$id <- 1

coefs <- do.call(rbind, lapply(list_m2[[2]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m2[[2]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids1)
newX <- model.matrix(list_m2[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pr1 <- as.data.frame(plogis(eta))
names(pr1)[1] <- "fit"
pr1$id <- 2

coefs <- do.call(rbind, lapply(list_m2[[3]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m2[[3]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids2)
newX <- model.matrix(list_m2[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pr2 <- as.data.frame(plogis(eta))
names(pr2)[1] <- "fit"
pr2$id <- 3

coefs <- do.call(rbind, lapply(list_m2[[4]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m2[[4]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] #else: returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids3)
newX <- model.matrix(list_m2[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pr3 <- as.data.frame(plogis(eta))
names(pr3)[1] <- "fit"
pr3$id <- 4

coefs <- do.call(rbind, lapply(list_m2[[5]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m2[[5]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids4)
newX <- model.matrix(list_m2[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pr4 <- as.data.frame(plogis(eta))
names(pr4)[1] <- "fit"
pr4$id <- 5

coefs <- do.call(rbind, lapply(list_m2[[6]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m2[[6]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids5)
newX <- model.matrix(list_m2[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pr5 <- as.data.frame(plogis(eta))
names(pr5)[1] <- "fit"
pr5$id <- 6

coefs <- do.call(rbind, lapply(list_m2[[7]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m2[[7]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
d1 <- as.vector(d1)
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
nd <- as.data.frame(ids6)
newX <- model.matrix(list_m2[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pr6 <- as.data.frame(plogis(eta))
names(pr6)[1] <- "fit"
pr6$id <- 7

coefs <- do.call(rbind, lapply(list_m2[[8]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m2[[8]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids7)
newX <- model.matrix(list_m2[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pr7 <- as.data.frame(plogis(eta))
names(pr7)[1] <- "fit"
pr7$id <- 8

coefs <- do.call(rbind, lapply(list_m2[[9]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m2[[9]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids8)
newX <- model.matrix(list_m2[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pr8 <- as.data.frame(plogis(eta))
names(pr8)[1] <- "fit"
pr8$id <- 9

coefs <- do.call(rbind, lapply(list_m2[[10]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m2[[10]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids9)
newX <- model.matrix(list_m2[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pr9 <- as.data.frame(plogis(eta))
names(pr9)[1] <- "fit"
pr9$id <- 10

list_yhat2 <- list(pr, pr1, pr2, pr3, pr4, pr5, pr6, pr7,pr8,pr9)

for (i in 1:10) {
f$yhat2[f$id==i] <- list_yhat2[[i]]$fit
}

#MODEL 3

coefs <- do.call(rbind, lapply(list_m3[[1]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m3[[1]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids)
newX <- model.matrix(list_m3[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pb <- as.data.frame(plogis(eta))
names(pb)[1] <- "fit"
pb$id <- 1

coefs <- do.call(rbind, lapply(list_m3[[2]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m3[[2]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids1)
newX <- model.matrix(list_m3[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pb1 <- as.data.frame(plogis(eta))
names(pb1)[1] <- "fit"
pb1$id <- 2

coefs <- do.call(rbind, lapply(list_m3[[3]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m3[[3]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids2)
newX <- model.matrix(list_m3[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pb2 <- as.data.frame(plogis(eta))
names(pb2)[1] <- "fit"
pb2$id <- 3

coefs <- do.call(rbind, lapply(list_m3[[4]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m3[[4]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] #else: returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids3)
newX <- model.matrix(list_m3[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pb3 <- as.data.frame(plogis(eta))
names(pb3)[1] <- "fit"
pb3$id <- 4

coefs <- do.call(rbind, lapply(list_m3[[5]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m3[[5]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids4)
newX <- model.matrix(list_m3[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pb4 <- as.data.frame(plogis(eta))
names(pb4)[1] <- "fit"
pb4$id <- 5

coefs <- do.call(rbind, lapply(list_m3[[6]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m3[[6]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids5)
newX <- model.matrix(list_m3[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pb5 <- as.data.frame(plogis(eta))
names(pb5)[1] <- "fit"
pb5$id <- 6

coefs <- do.call(rbind, lapply(list_m3[[7]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m3[[7]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
d1 <- as.vector(d1)
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
nd <- as.data.frame(ids6)
newX <- model.matrix(list_m3[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pb6 <- as.data.frame(plogis(eta))
names(pb6)[1] <- "fit"
pb6$id <- 7

coefs <- do.call(rbind, lapply(list_m3[[8]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m3[[8]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids7)
newX <- model.matrix(list_m3[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pb7 <- as.data.frame(plogis(eta))
names(pb7)[1] <- "fit"
pb7$id <- 8

coefs <- do.call(rbind, lapply(list_m3[[9]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m3[[9]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids8)
newX <- model.matrix(list_m3[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pb8 <- as.data.frame(plogis(eta))
names(pb8)[1] <- "fit"
pb8$id <- 9

coefs <- do.call(rbind, lapply(list_m3[[10]], function(i) coef(summary(i))[,1]))
ses <- do.call(rbind, lapply(list_m3[[10]], function(i) coef(summary(i))[,2]))
d1 <- as.matrix(t(do.call(rbind, mi.meld(coefs, ses))))
d1<- d1[,-2]
#d1 <- d1[-3] use if returns "Error in newX %*% d1 : non-conformable arguments"
d1 <- as.vector(d1)
nd <- as.data.frame(ids9)
newX <- model.matrix(list_m3[[1]][[1]]$fixed,nd)
eta <- newX %*% d1
pb9 <- as.data.frame(plogis(eta))
names(pb9)[1] <- "fit"
pb9$id <- 10

# merge predictions 

list_yhat3 <- list(pb, pb1, pb2, pb3, pb4, pb5, pb6, pb7,pb8,pb9)
 

# integrate predictions into the dataset

for (i in 1:10) {
f$yhat3[f$id==i] <- list_yhat3[[i]]$fit
}

# ROC PLOT OR FIGURE 4 IN THE ARTICLE

fd <- f

pred <- prediction(abs(fd$yhat), fd$ally)
perf <- performance(pred, "tpr", "fpr")
auc <- performance(pred,"auc")
# now converting S4 class to vector
auc <- unlist(slot(auc, "y.values"))

pred2 <- prediction(abs(fd$yhat2), fd$ally)
perf2 <- performance(pred2, "tpr", "fpr")
auc2 <- performance(pred2,"auc")
# now converting S4 class to vector
auc2 <- unlist(slot(auc2, "y.values"))

pred3 <- prediction(abs(fd$yhat3), fd$ally)
perf3 <- performance(pred3, "tpr", "fpr")
auc3 <- performance(pred3,"auc")
# now converting S4 class to vector
auc3 <- unlist(slot(auc3, "y.values"))

pdf('ROC.pdf', 7, 7)
  par(xpd=F, mar=c(5,5,1,3), mfrow=c(2,2))
par(mar=c(5,5,2,2), xaxs = "i",yaxs = "i",cex.axis=1.3,cex.lab=1.4)
plot(perf,col='grey20',main='Model 1')
text(.75,.4,paste('AUC = ',round(auc,2)),col='grey20')
plot(perf2,col='grey40',main='Model 2')
text(.75,.4,paste('AUC = ',round(auc2,2)),col='grey40')
plot(perf3,col='grey60',main='Model 3')
text(.75,.4,paste('AUC = ',round(auc3,2)),col='grey60')
dev.off()

# Model Criticism Plots or FIGURE 5 IN THE ARTICLE
# add dyad and year

fd <- cbind(d$dyad, d$year, fd)
names(fd)[1] <- "dyad"
names(fd)[2] <- "year"
names(fd)[3] <- "ally"

# model 1 (shown in the article as FIGURE 5)
m1diag<-DiagPlot(
	f=fd$yhat,
	y=fd$ally,
	labels= paste(fd$dyad,fd$year,sep='-'),
	worstN=10,
	label_spacing = 25,
	lab_adjust=.4,
	right_margin=11,
	top_margin=5,
	text_size=6,
	title=""
	)
	
# model 2 (FIGURE 4 IN THE APPENDIX)
m2diag<-DiagPlot(
	f=fd$yhat2,
	y=fd$ally,
	labels= paste(fd$dyad,fd$year,sep='-'),
	worstN=10,
	label_spacing = 25,
	lab_adjust=.4,
	right_margin=11,
	top_margin=4,
	text_size=6,
	title=""
	)

# model 3 (FIGURE 5 IN THE APPENDIX)
m3diag<-DiagPlot(
	f=fd$yhat3,
	y=fd$ally,
	labels= paste(fd$dyad,fd$year,sep='-'),
	worstN=10,
	label_spacing = 25,
	lab_adjust=.4,
	right_margin=11,
	top_margin=4,
	text_size=6,
	title=""
	)
	
pdf("criticm1.pdf")
grid.draw(m1diag)
dev.off()

pdf("criticm2.pdf")
grid.draw(m2diag)
dev.off()

pdf("criticm3.pdf")
grid.draw(m3diag)
dev.off()

#################################################
### APPENDIX ####################################
#################################################

# mean predicted probabilities for each dyad based on imputed data (FIGURE 1 IN THE APPENDIX)

p_eta <- lapply(g, predict, newdata=b$Imputations, re.form=NA, type = "link", se = TRUE)
pr_eta <- do.call(rbind, p_eta)
l = apply( cbind( pr_eta ) , 1 , unlist )
fdd <- t(l[c(1:985),]) 
sdd <- t(l[c(986:1970),] )
da <- data.frame(t(do.call(rbind, mi.meld(fdd, sdd))))
names(da)[1] <- "prob"
names(da)[2] <- "se"
da$left <- plogis(da$prob-(0.67449*da$se))
da$right <- plogis(da$prob+(0.67449*da$se))
da$fit <- plogis(da$prob)
dfd <- data.frame(da, d)

ag <- ddply(dfd, .(dyadid, dyad), summarise, prob = mean(fit), lo=mean(left), up=mean(right), sponsor=max(sponsor))
ag$dyad <- toupper(ag$dyad)
ag$sponsor <- as.factor(ag$sponsor)

levels(ag$sponsor) <- gsub("0", "No Sponsor", levels(ag$sponsor))
levels(ag$sponsor) <- gsub("1", "Single Sponsor", levels(ag$sponsor))
levels(ag$sponsor) <- gsub("2", "Different Sponsor", levels(ag$sponsor))
levels(ag$sponsor) <- gsub("3", "Shared Sponsor", levels(ag$sponsor))

ag$color[ag$sponsor=="No Sponsor"] <- "grey80"
ag$color[ag$sponsor=="Different Sponsor"] <- "grey60"
ag$color[ag$sponsor=="Single Sponsor"] <- "grey40"
ag$color[ag$sponsor=="Shared Sponsor"] <- "grey20"

ag$dyad2 <- reorder(ag$dyad, ag$prob)

mp <- ggplot(data = ag, aes(y =prob, x=interaction(dyad2,sponsor), ymin = lo, ymax = up, colour = factor(sponsor),  shape=factor(sponsor))) +
  scale_x_discrete("dyad",breaks=interaction(ag$dyad,ag$sponsor),labels=ag$dyad) +
  scale_y_continuous(breaks=c(0,.2,.4,.6,.8,1)) +
  geom_point(position = position_dodge(width = 0.2), size=3) +
  geom_errorbar(position = position_dodge(width = 0.2), width = 0.1, size=1.5) +
  geom_text(aes(y =-2.5, x=interaction(ag$dyad2,ag$sponsor), label=ag$dyad, hjust=0, color=ag$sponsor), size=5) +
  coord_flip() +
  scale_colour_manual(values = c("grey80", "grey60", "grey40", "grey20")) +
  scale_shape_manual(values = c(16,15,17,16)) +
  theme_bw() + xlab("") + ylab("") +
  theme(panel.background = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major=element_blank(),
        panel.border=element_blank(), 
        axis.ticks=element_blank(), 
        axis.text.y = element_blank(),
        axis.text.x = element_text(size=14),	
        axis.line.y=element_blank(),
        axis.line.x=element_blank(),
        axis.title.x=element_blank(), 
        axis.title.y = element_blank(),
        legend.title = element_blank(),
        legend.key=element_rect(fill=NA),
        legend.position="right", 
        legend.direction="vertical",
		legend.text=element_text(size=16)) + 
        guides(colour = guide_legend(override.aes = list(size=5, linetype=0, shape=c(16,15,17,16))))

ggsave("fig5.pdf", height=50*(gr^-1), width=12.5, mp)

# confusion matrix (FIGURE 2 IN THE APPENDIX)

fd$ally <- as.factor(fd$ally)

obs <- as.vector(fd$ally)
predicted <- as.vector(fd$yhat)
confusion.matrix(obs, predicted, threshold = 0.5)
accuracy(obs, predicted, threshold = 0.5)

obs2 <- as.vector(fd$ally)
predicted2 <- as.vector(fd$yhat2)
confusion.matrix(obs2, predicted2, threshold = 0.5)
accuracy(obs2, predicted2, threshold = 0.5)

obs3 <- as.vector(fd$ally)
predicted3 <- as.vector(fd$yhat3)
confusion.matrix(obs3, predicted3, threshold = 0.5)
accuracy(obs3, predicted3, threshold = 0.5)

cm1 <- matrix(c(700,  91, 79,  115),ncol=2,byrow=TRUE)
colnames(cm1) <- c("Actual = 0","Actual = 1")
rownames(cm1) <- c("Predicted = 0","Predicted = 1")

cm2 <- matrix(c(701,  94, 78,  112),ncol=2,byrow=TRUE)
colnames(cm2) <- c("Actual = 0","Actual = 1")
rownames(cm2) <- c("Predicted = 0","Predicted = 1")

cm3 <- matrix(c(703,  86, 76,  120),ncol=2,byrow=TRUE)
colnames(cm3) <- c("Actual = 0","Actual = 1")
rownames(cm3) <- c("Predicted = 0","Predicted = 1")

pdf('ConfusionMatrix.pdf', 8, 7)
       par(xpd=F, mar=c(5,5,1,3), mfrow=c(2,2))	
             fourfoldplot(cm1,
             main = "",
             color=c("grey20", "grey80"))
			 title(main="Model 1\n Proportion correct = 0.82", cex.main=1.05, line=-2)
			 fourfoldplot(cm2,
             main = "",
             color=c("grey20", "grey80"))
			 title(main="Model 2\n Proportion correct = 0.82", cex.main=1.05, line=-2)
			 fourfoldplot(cm3,
             main = "",
             color=c("grey20", "grey80"))
			 title(main="Model 3\n Proportion correct = 0.83", cex.main=1.05, line=-2)
dev.off()	

## Precision-recall (PR) plot or FIGURE 3 IN THE APPENDIX

fp1 <- prediction(fd$yhat,fd$ally)
fplot1 <- performance(fp1,'prec','rec')
fp2 <- prediction(fd$yhat2,fd$ally)
fplot2 <- performance(fp2,'prec','rec')
fp3 <- prediction(fd$yhat3,fd$ally)
fplot3 <- performance(fp3,'prec','rec')

pdf('precision.pdf')
par(mfrow=c(2,2))
plot(fplot1,col='firebrick2',main='Model 1',ylim=c(0,1))
plot(fplot2,col='dodgerblue2',main='Model 2',ylim=c(0,1))
plot(fplot3,col='green4',main='Model 3',ylim=c(0,1))

plot(fplot1,col='firebrick2',ylim=c(0,1),main='Joined')
plot(fplot2,col='dodgerblue2',add=T,ylim=c(0,1))
plot(fplot3,col='green4',add=T,ylim=c(0,1))
par(mfrow=c(1,1))
dev.off()

#### ROBUSTNESS --- POST-2001 FIGURE 6 IN THE APPENDIX####

ad <- a
var_out <- c("insurgentsa", "insurgentsb", 'size_a',	'size_b',	'id',	'location',	'oldid',	'target', 'dur')

fdd <- ad[, !(names(ad)%in%var_out)]
pw = fdd[!fdd$year > 2001,]

pw$cowcode <- as.factor(pw$cowcode)

pw$ratio[pw$ratio==1] <- 1-1e-3
pw$ratio[pw$ratio==0] <- 0+1e-3

pw$relfrac[pw$relfrac==1] <- 1-1e-3
pw$relfrac[pw$relfrac==0] <- 0+1e-3

pw$ef[pw$ef==1] <- 1-1e-3
pw$ef[pw$ef==0] <- 0+1e-3

set.seed(1234)

wp <- amelia(x=pw, m=5, 
ts='year', 
cs='cowcode', 
lgstc='ratio', 
noms='wd', 
idvars=c('dyad', 'dyadid'),
bounds=rbind(c(10, 0, 1), c(9, 0, 1), c(11, 0, Inf), c(12, 0, Inf), c(17, 0, Inf)),
polytime=3)

f0 <- ally~ 1+statesupport*wd+lgdppc+lnmilex+ldur+ncontig+relfrac+ef+lndurable+(1|cowcode)

f1 <- ally~ 1+factor(sponsor)+lgdppc+lnmilex+ldur+ncontig+relfrac+ef+lndurable+(1|cowcode)

f0b <- ally~ 1+rsup*wd+lgdppc+lnmilex+ldur+ncontig+relfrac+ef+lndurable+(1|cowcode)

n0 <- lapply(wp$imputations, burrito, form=f0)
n1 <- lapply(wp$imputations, burrito, form=f1)
m0 <- lapply(wp$imputations, burrito, form=f0b)

#save(m0, file = "post2001model2.rda")
#save(n0, file = "post2001model2bbsupport.rda")
#save(g0b, file = "post2001model1.rda")

save(list=c('m0', 'n0', 'n1'), file='post2001models2.RData')

# pool coefficients for model 1 by Rubin's rule

gcoefs6 <- do.call(rbind, lapply(n1, function(i) coef(summary(i))[,1])) # coef estimates
gses6 <- do.call(rbind, lapply(n1, function(i) coef(summary(i))[,2])) # standard errors
gd6 <- data.frame(t(do.call(rbind, mi.meld(gcoefs6, gses6)))) # pool coefs and errors by Rubin's rule 
names(gd6) <- c("coef", "se")
gd6 <- gd6[-c(1), ] # get rid of intercept
gd6$lo <- gd6$coef - 1.96*gd6$se # lower 95% confidence intervals
gd6$up <- gd6$coef + 1.96*gd6$se # upper 95% confidence intervals
gd6$var = c("Single sponsor", "Different sponsor", "Shared sponsor", "GDP p.c. (ln)", "Expenditure (ln)", "Duration (ln)", "Non-Contiguity", "Religious frac.", "Ethnic frac.", "Durability (ln)") # assign names to vars

# plot model 7

p6 <- ggplot(data=gd6) + scale_x_discrete(limits = c("Durability (ln)", "Ethnic frac.", "Religious frac.", "Non-Contiguity", "Duration (ln)", "Expenditure (ln)", "GDP p.c. (ln)", "Shared sponsor", "Different sponsor",   "Single sponsor"))
p6 <- p6 + geom_hline(yintercept = 0, colour = "grey20", lty = 2) 
p6 <- p6 + geom_linerange(aes(x = var, ymin = lo,
                            ymax = up), col= "grey20",
                            lwd = 1, position = position_dodge(width = 1/2))
p6 <- p6 + geom_pointrange(aes(x = var, y = coef, ymin = lo,
                                 ymax = up), col= "grey20",
                             lwd = 1/2, position = position_dodge(width = 1/2),
                             shape = 21, fill = "grey20")
p6 <- p6 + coord_flip() + guides(fill=FALSE) + theme(panel.background = element_blank(),
              panel.grid.minor = element_blank(),
              panel.grid.major=element_blank(),
              panel.border=element_blank(),
			  plot.title=element_text(hjust = .5, size=11),
              axis.title.x = element_text(colour="grey30"),
              axis.title.y = element_text(colour="grey30"),
              axis.line   = element_line(colour=NA),
              axis.line.x = element_line(colour="grey30"),
              axis.line.y = element_line(colour="grey30"),
              axis.text.x = element_text(colour="grey30", hjust = 1),
              axis.text.y = element_text(colour="grey30"), legend.title = element_blank(),
              legend.key=element_rect(fill=NA),
              legend.position="bottom") + ylab("Coefficient Estimates") + xlab(NULL) +
			  ggtitle("Model 7")

# pool coefficients for model 8 by Rubin's rule

gcoefs7 <- do.call(rbind, lapply(m0, function(i) coef(summary(i))[,1]))
gses7 <- do.call(rbind, lapply(m0, function(i) coef(summary(i))[,2]))
gd7 <- data.frame(t(do.call(rbind, mi.meld(gcoefs7, gses7))))
names(gd7) <- c("coef", "se") #change names
gd7 <- gd7[-c(1), ] # get rid of intercept
gd7$lo <- gd7$coef - 1.96*gd7$se # lower 95% confidence intervals
gd7$up <- gd7$coef + 1.96*gd7$se # upper 95% confidence intervals
gd7$var = c("Sponsor", "Weak Dyad", "GDP p.c. (ln)", "Expenditure (ln)", "Duration (ln)", "Non-Contiguity", "Religious frac.", "Ethnic frac.", "Durability (ln)", "Sponsor x\n Weak Dyad")

# plot model 2

p7 <- ggplot(data=gd7) + scale_x_discrete(limits = c("Durability (ln)", "Ethnic frac.", "Religious frac.", "Non-Contiguity", "Duration (ln)", "Expenditure (ln)", "GDP p.c. (ln)","Sponsor x\n Weak Dyad", "Sponsor", "Weak Dyad"))
p7 <- p7 + geom_hline(yintercept = 0, colour = "grey20", lty = 2) 
p7 <- p7 + geom_linerange(aes(x = var, ymin = lo,
                            ymax = up), col= "grey20",
                            lwd = 1, position = position_dodge(width = 1/2))
p7 <- p7 + geom_pointrange(aes(x = var, y = coef, ymin = lo,
                                 ymax = up), col= "grey20",
                             lwd = 1/2, position = position_dodge(width = 1/2),
                             shape = 21, fill = "grey20")
p7 <- p7 + coord_flip() + guides(fill=FALSE) + theme(panel.background = element_blank(),
              panel.grid.minor = element_blank(),
              panel.grid.major=element_blank(),
              panel.border=element_blank(),
			  plot.title=element_text(hjust = .5, size=11),
              axis.title.x = element_text(colour="grey30"),
              axis.title.y = element_text(colour="grey30"),
              axis.line   = element_line(colour=NA),
              axis.line.x = element_line(colour="grey30"),
              axis.line.y = element_line(colour="grey30"),
              axis.text.x = element_text(colour="grey30", hjust = 1),
              axis.text.y = element_text(colour="grey30"), legend.title = element_blank(),
              legend.key=element_rect(fill=NA),
              legend.position="bottom") + ylab("Coefficient Estimates") + xlab(NULL) +
			  ggtitle("Model 8")

# COMBINE MODELS INTO A SINGLE PLOT
			  
pp7 <- grid.arrange(p6, p7, ncol=2, nrow=1)			  
ggsave("fig7.pdf", height=7, width=7, pp7)

# plot model 9 (FIGURE 7)

gcoefs8 <- do.call(rbind, lapply(n0, function(i) coef(summary(i))[,1]))
gses8 <- do.call(rbind, lapply(n0, function(i) coef(summary(i))[,2]))
gd8 <- data.frame(t(do.call(rbind, mi.meld(gcoefs8, gses8))))
names(gd8) <- c("coef", "se") #change names
gd8 <- gd8[-c(1), ] # get rid of intercept
gd8$lo <- gd8$coef - 1.96*gd8$se # lower 95% confidence intervals
gd8$up <- gd8$coef + 1.96*gd8$se # upper 95% confidence intervals
gd8$var = c("Sponsor", "Weak Dyad", "GDP p.c. (ln)", "Expenditure (ln)", "Duration (ln)", "Non-Contiguity", "Religious frac.", "Ethnic frac.", "Durability (ln)", "Sponsor x\n Weak Dyad")

# plot model 2

p8 <- ggplot(data=gd8) + scale_x_discrete(limits = c("Durability (ln)", "Ethnic frac.", "Religious frac.", "Non-Contiguity", "Duration (ln)", "Expenditure (ln)", "GDP p.c. (ln)","Sponsor x\n Weak Dyad", "Sponsor", "Weak Dyad"))
p8 <- p8 + geom_hline(yintercept = 0, colour = "grey20", lty = 2) 
p8 <- p8 + geom_linerange(aes(x = var, ymin = lo,
                            ymax = up), col= "grey20",
                            lwd = 1, position = position_dodge(width = 1/2))
p8 <- p8 + geom_pointrange(aes(x = var, y = coef, ymin = lo,
                                 ymax = up), col= "grey20",
                             lwd = 1/2, position = position_dodge(width = 1/2),
                             shape = 21, fill = "grey20")
p8 <- p8 + coord_flip() + guides(fill=FALSE) + theme(panel.background = element_blank(),
              panel.grid.minor = element_blank(),
              panel.grid.major=element_blank(),
              panel.border=element_blank(),
			  plot.title=element_text(hjust = .5, size=11),
              axis.title.x = element_text(colour="grey30"),
              axis.title.y = element_text(colour="grey30"),
              axis.line   = element_line(colour=NA),
              axis.line.x = element_line(colour="grey30"),
              axis.line.y = element_line(colour="grey30"),
              axis.text.x = element_text(colour="grey30", hjust = 1),
              axis.text.y = element_text(colour="grey30"), legend.title = element_blank(),
              legend.key=element_rect(fill=NA),
              legend.position="bottom") + ylab("Coefficient Estimates") + xlab(NULL) +
			  ggtitle("Model 9")

ggsave("fig8.pdf", height=7, width=7, p8)

# Cold War Period for rebel-level variables (FIGURE 8 IN THE APPENDIX))

# Create Cold War var

d$cw[d$year<=1990] <- 1
d$cw[d$year>=1991] <- 0

# run a model without interaction

r3 <- glmmadmb(ally~1+ideol+splinter+ethnic+ldur+cw+(1|cowcode), family = "binomial", data=d)

# run a model with interaction between shared ideology and period

r4 <- glmmadmb(ally~1+ideol*cw+splinter+ethnic+ldur+(1|cowcode), family = "binomial", data=d)

# make plot

dr3 <- data.frame(Variable = rownames(summary(r3)$coef),
                          Coefficient = summary(r3)$coef[, 1],
                          SE = summary(r3)$coef[, 2])
dr4 <- data.frame(Variable = rownames(summary(r4)$coef),
                          Coefficient = summary(r4)$coef[, 1],
                          SE = summary(r4)$coef[, 2])
						  
dr3$modelName <- "Model 10: Cold War"
dr4$modelName <- "Model 11: Interaction"

dr <- data.frame(rbind(dr3, dr4))

# change name of variables for the plot

levels(dr$Variable) <- gsub("ideol", "Shared ideology", levels(dr$Variable))
levels(dr$Variable) <- gsub("splinter", "Splinter", levels(dr$Variable))
levels(dr$Variable) <- gsub("ethnic", "Shared ethnic pool", levels(dr$Variable))
levels(dr$Variable) <- gsub("ldur", "Duration (ln)", levels(dr$Variable))
levels(dr$Variable) <- gsub("cw", "Cold War", levels(dr$Variable))
levels(dr$Variable) <- gsub("ideol:cw", "Ideology x\n Cold War", levels(dr$Variable))

# discard intercepts

dr <- dr[-c(1, 7), ]

ci.95 <- -qnorm((1-0.95)/2)  # 95% multiplier


# make plot

posr3 <- c( "Duration (ln)", "Shared ethnic pool", "Splinter", "Shared ideology:Cold War", "Cold War","Shared ideology")


pd <- ggplot(data=dr, aes(shape = modelName, col=modelName)) + scale_x_discrete(limits = posr3)
pd <- pd + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
pd <- pd + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*ci.95,
                                ymax = Coefficient + SE*ci.95),
                            lwd = 1, position = position_dodge(width = 0.95))
pd <- pd + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*ci.95,                    
ymax = Coefficient + SE*ci.95),
lwd = 0.55, position = position_dodge(width = 0.95)) + 
scale_shape_manual(values=c(16,17)) + 
scale_color_manual(values=c("grey20", "grey80"))
pd <- pd + coord_flip() + guides(colour = guide_legend(override.aes = list(size=1, linetype=0, shape=c(16,17)))) +
theme(panel.background = element_blank(),
              panel.grid.minor = element_blank(),
              panel.grid.major=element_blank(),
              panel.border=element_blank(),
              axis.title.x = element_text(colour="grey30"),
              axis.title.y = element_text(colour="grey30"),
              axis.line   = element_line(colour=NA),
              axis.line.x = element_line(colour="grey30"),
              axis.line.y = element_line(colour="grey30"),
              axis.text.x = element_text(colour="grey30", hjust = 1),
              axis.text.y = element_text(colour="grey30"), legend.title = element_blank(),
              legend.key=element_rect(fill=NA),
              legend.position="bottom", legend.direction="vertical") + 
              xlab(NULL) + 
              ylab("Coefficient Estimates")

 ggsave("period.pdf", width=7, height=7, pd)

######################
	### Make tables with 
	### standard errors 
	### and p-values TABLE 1 in the Appendix
	######################
# pool AIC and BIC

# base

aic0 <- lapply(base, extractAIC)
bic0 <- lapply(base, BIC)
ai0 <- do.call(cbind, aic0)
bi0 <- do.call(cbind, bic0)
l0 = apply( rbind( bi0 ) , 1 , unlist )
l0b = apply( rbind( ai0 ) , 1 , unlist )
aic0 <- as.matrix(l0b[,c(2)])
bic0 <- as.matrix(l0[,c(1)])
adf0 <- data.frame(t(do.call(rbind, mi.meld(bic0, aic0))))

# model 1

aic1 <- lapply(g, extractAIC)
bic1 <- lapply(g, BIC)
ai1 <- do.call(cbind, aic1)
bi1 <- do.call(cbind, bic1)
l1 = apply( rbind( bi1 ) , 1 , unlist )
l1b = apply( rbind( ai1 ) , 1 , unlist )
aic1 <- as.matrix(l1b[,c(2)])
bic1 <- as.matrix(l1[,c(1)])
adf1 <- data.frame(t(do.call(rbind, mi.meld(bic1, aic1))))

# model 2

aic2<- lapply(g0, extractAIC)
bic2 <- lapply(g0, BIC)
ai2 <- do.call(cbind, aic2)
bi2 <- do.call(cbind, bic2)
l2 = apply( rbind( bi2 ) , 1 , unlist )
l2b = apply( rbind( ai2 ) , 1 , unlist )
aic2 <- as.matrix(l2b[,c(2)])
bic2 <- as.matrix(l2[,c(1)])
adf2 <- data.frame(t(do.call(rbind, mi.meld(bic2, aic2))))

# model 3

aic3 <- lapply(g0b, extractAIC)
bic3 <- lapply(g0b, BIC)
ai3 <- do.call(cbind, aic3)
bi3 <- do.call(cbind, bic3)
l3 = apply( rbind( bi3 ) , 1 , unlist )
l3b = apply( rbind( ai3 ) , 1 , unlist )
aic3 <- as.matrix(l3b[,c(2)])
bic3 <- as.matrix(l3[,c(1)])
adf3 <- data.frame(t(do.call(rbind, mi.meld(bic3, aic3))))

# pool Log-Likelihood
# Average Log-Likelihood

llbase <- do.call(rbind, lapply(base, function(i) logLik(i)))
colMeans(llbase)
llg0b <- do.call(rbind, lapply(g0b, function(i) logLik(i)))
colMeans(llg0b)
llg0 <- do.call(rbind, lapply(g0, function(i) logLik(i)))
colMeans(llg0)
llg <- do.call(rbind, lapply(g, function(i) logLik(i)))
colMeans(llg)

# pool coefficients, standard errors and calculate p-values
	
fits0 <- MIcombine(base)
md0 <- summary(fits0)
mdd0 <- cbind(md0, df = fits0$df)
mdd0$z <- mdd0$results/mdd0$se
mdd0$pvalue <- exp(-0.717*mdd0$z - 0.416*mdd0$z^2)	
	
fits1 <- MIcombine(g)
md <- summary(fits1)
mdd <- cbind(md, df = fits1$df)
mdd$z <- mdd$results/mdd$se
mdd$pvalue <- exp(-0.717*mdd$z - 0.416*mdd$z^2)

fits2 <- MIcombine(g0)
md2 <- summary(fits2)
mdd2 <- cbind(md2, df = fits2$df)
mdd2$z <- mdd2$results/mdd2$se
mdd2$pvalue <- exp(-0.717*mdd2$z - 0.416*mdd2$z^2)

fits3 <- MIcombine(g0b)
md3 <- summary(fits3)
mdd3 <- cbind(md3, df = fits3$df)
mdd3$z <- mdd3$results/mdd3$se
mdd3$pvalue <- exp(-0.717*mdd3$z - 0.416*mdd3$z^2)

### Robustness 2 - use rebel level variables as controls (TABLE 2 IN THE APPENDIX)

d$dyadid <- as.factor(d$dyadid)

r1 <- glmmadmb(ally~1+ideol+splinter+ethnic+ldur+(1|cowcode), family = "binomial", data=d)
r2 <- glmmadmb(ally~1+factor(sponsor)+ideol+splinter+ethnic+ldur+(1|cowcode), family = "binomial", data=d)

summary(r1)
summary(r2)

# Heckman model (TABLE 3 IN THE APPENDIX)

attach(d)

selection1 <- selection(selection = rsup ~ 1 + lgdppc + lnmilex + wd + relfrac + ef, outcome = ally ~ lndurable + ldur + ideol + ethnic + splinter, 
data = d, method = "2step")
summary(selection1)

##############################################################################################################################################
#### END REPLICATION
##############################################################################################################################################